# Import data from .RData
#load("data/PSZ_cont_snotel.RData")
load("data/cont_snotel.RData")
# Count the number of SNOTEL Sites per Ecoregion
eco_sntl_sites <- cont_snotel %>% 
  group_by(US_L3NAME, US_L3CODE) %>% 
  summarize(eco_sntl_sites = n_distinct(site_id)) %>% 
  filter(eco_sntl_sites >= 5)
## `summarise()` has grouped output by 'US_L3NAME'. You can override using the
## `.groups` argument.
# All SNOTEL Sites, Calculate spring days with new SWE and new SWE depth 
# filter by date (March 1st - 60/61, March 15th - 74/75, April 1st - 91/92, April 15th - 105/106, August 1st - 243/244)
# Filter data by date from 2000 to present and from March 1st to August 1st
# group by year and site
start_year = 2000
# start_day = 60 #March 1
# cont_snotel_spring_mar1 <- cont_snotel %>%
#   mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
#   mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
#   filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243, 
#                                                    yday >= start_day+1 & yday <= 243+1)) %>% 
#   select(-c("X", "network", "description", "start", "end"))
# 
# start_day = 74 #March 15
# cont_snotel_spring_mar15 <- cont_snotel %>%
#   mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
#   mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
#   filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243, 
#                                                    yday >= start_day+1 & yday <= 243+1)) %>% 
#   select(-c("X", "network", "description", "start", "end"))

start_day = 91 #April 1
cont_snotel_spring_apr1 <- cont_snotel %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
  mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent),
                         snow_water_equivalent-lag(snow_water_equivalent), 0),
         newSWE = ifelse(temperature_min > 0, 0, newSWE)) %>%
  filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243, 
                                                   yday >= start_day+1 & yday <= 243+1)) %>% 
  select(-c("X", "network", "description", "start", "end"))

# start_day = 105 #April 15
# cont_snotel_spring_apr15 <- cont_snotel %>%
#   mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
#   mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent)+3, snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
#   filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243, 
#                                                    yday >= start_day+1 & yday <= 243+1)) %>% 
#   select(-c("X", "network", "description", "start", "end"))
# # Determine number of days with increased SWE per spring at each site (YEARLY - March 1)
# cont_snotel_site_yr_mar1 <- cont_snotel_spring_mar1 %>%
#   group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
#   summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>% 
#   merge(., eco_sntl_sites, by = c("US_L3NAME")) %>% 
#   mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "3/1")
# 
# # Determine number of days with increased SWE per spring at each site (YEARLY - March 15)
# cont_snotel_site_yr_mar15 <- cont_snotel_spring_mar15 %>%
#   group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
#   summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>% 
#   merge(., eco_sntl_sites, by = c("US_L3NAME")) %>% 
#   mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "3/15")

# Determine number of days with increased SWE per spring at each site (YEARLY - April 1)
cont_snotel_site_yr_apr1 <- cont_snotel_spring_apr1 %>%
  group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
  summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>% 
  merge(., eco_sntl_sites, by = c("US_L3NAME", "US_L3CODE")) %>% 
  mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "4/1") 
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude', 'elev', 'US_L3CODE', 'US_L3NAME'. You can override
## using the `.groups` argument.
eco_sntl_sites_med <- cont_snotel_site_yr_apr1 %>% 
  group_by(eco_sntl) %>% 
  summarize(median_spring_SWE = median(newSWE_days, na.rm = T))

# # Determine number of days with increased SWE per spring at each site (YEARLY - April 15)
# cont_snotel_site_yr_apr15 <- cont_snotel_spring_apr15 %>%
#   group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
#   summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>% 
#   merge(., eco_sntl_sites, by = c("US_L3NAME")) %>% 
#   mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "4/15")

# Merge the DF's into one large datafame for plotting
# cont_snotel_site_yr <- rbind(cont_snotel_site_yr_mar1, cont_snotel_site_yr_mar15, cont_snotel_site_yr_apr1, cont_snotel_site_yr_apr15)

# Determine average number of days with increased SWE per spring at each site in the PSZ (SITE - April 1)
cont_snotel_site_apr1 <- cont_snotel_site_yr_apr1 %>%
  group_by(site_id, site_name, state, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY, eco_sntl_sites, eco_sntl) %>%
  summarise(mean_spring_days = mean(newSWE_days), med_spring_days = median(newSWE_days),
            mean_newSWE = mean(newSWE), med_newSWE = median(newSWE)) %>% 
  merge(.,eco_sntl_sites_med, by = "eco_sntl")
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude', 'elev', 'US_L3CODE', 'US_L3NAME', 'L3_KEY',
## 'eco_sntl_sites'. You can override using the `.groups` argument.
# # View data
# mapview(PSZ_cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "med_spring_days",
#         layer.name = "AVG Sping Days (PSZ Sites)" , crs = 4326, grid = FALSE) +
#   mapview(eco_L3_4326, zcol = "L3_KEY")
# Spring SWE increases (After April 1)
spring_snow_days <- ggplot(cont_snotel_site_yr_apr1, aes(x = reorder(US_L3CODE,newSWE_days,FUN = median, na.rm = TRUE), y= newSWE_days, 
                                                         fill = reorder(eco_sntl,newSWE_days,FUN = median,na.rm = TRUE))) +
  geom_boxplot() +
  labs(x = "EPA Level III Ecoregion", y = "Days of Increased SWE", fill = "") +
  scale_y_continuous(breaks = seq(0,50,10), limits = c(0,50), expand = c(0.01, 0.01)) +
  scale_x_discrete(labels = function(eco_sntl) str_wrap(eco_sntl, width = 28)) +
  theme_bw()
spring_snow_days

ggplotly(spring_snow_days)
ggplot_build(spring_snow_days)$data
## [[1]]
##       fill ymin lower middle upper ymax
## 1  #F8766D    0     0      0     1    2
## 2  #E58700    0     1      4     9   21
## 3  #C99800    0     0      4     9   22
## 4  #A3A500    0     1      4     9   21
## 5  #6BB100    0     1      5    10   23
## 6  #00BA38    0     3      6    10   20
## 7  #00BF7D    0     2      6    12   27
## 8  #00C0AF    0     3      7    13   28
## 9  #00BCD8    0     3      7    12   23
## 10 #00B0F6    0     2      7    12   27
## 11 #619CFF    0     3      7    13   28
## 12 #B983FF    0     3      8    13   28
## 13 #E76BF3    0     5     10    16   32
## 14 #FD61D1    0     6     11    18   36
## 15 #FF67A4    0     8     14    20   38
##                                                                                                                                                               outliers
## 1  3, 6, 3, 4, 5, 6, 7, 3, 6, 3, 5, 5, 6, 13, 4, 3, 3, 6, 6, 3, 3, 5, 4, 11, 4, 4, 7, 3, 6, 4, 10, 4, 6, 4, 3, 3, 7, 3, 3, 6, 3, 3, 5, 6, 5, 3, 3, 3, 6, 4, 3, 3, 3, 4
## 2                                           23, 24, 25, 29, 37, 32, 24, 27, 29, 26, 25, 38, 25, 25, 33, 34, 27, 22, 30, 24, 33, 29, 27, 25, 24, 38, 24, 22, 32, 26, 23
## 3                                                                                                                               24, 23, 25, 23, 25, 24, 24, 23, 24, 32
## 4                                                                                                                                   27, 23, 24, 23, 26, 27, 22, 28, 26
## 5                                                                                                                               28, 24, 30, 26, 37, 33, 24, 27, 25, 28
## 6                                                                                               27, 23, 25, 25, 28, 27, 21, 27, 25, 33, 29, 27, 22, 23, 35, 22, 21, 21
## 7                                                                                                       31, 36, 30, 29, 28, 35, 29, 35, 28, 30, 28, 32, 35, 35, 29, 33
## 8                                                                                                                                                       29, 36, 32, 33
## 9                                                                                                                                                                     
## 10                                                                                          29, 31, 30, 36, 29, 33, 31, 28, 30, 30, 34, 28, 29, 29, 28, 30, 31, 29, 29
## 11                                  33, 31, 34, 38, 33, 29, 40, 40, 30, 34, 33, 29, 30, 30, 31, 31, 29, 29, 43, 30, 32, 30, 30, 33, 30, 34, 30, 29, 29, 29, 30, 30, 34
## 12                                                                                                                                                  29, 31, 29, 32, 29
## 13                                                                                      34, 39, 36, 37, 33, 34, 35, 35, 37, 36, 34, 35, 34, 38, 37, 34, 34, 38, 43, 34
## 14                                                                                                                                                                  38
## 15                                                                                                                      39, 39, 41, 39, 39, 40, 41, 43, 41, 41, 48, 41
##     notchupper  notchlower  x flipped_aes PANEL group ymin_final ymax_final
## 1   0.06366001 -0.06366001  1       FALSE     1     1          0         13
## 2   4.48401062  3.51598938  2       FALSE     1     2          0         38
## 3   4.55818383  3.44181617  3       FALSE     1     3          0         32
## 4   4.58806602  3.41193398  4       FALSE     1     4          0         28
## 5   5.53593641  4.46406359  5       FALSE     1     5          0         37
## 6   6.43214996  5.56785004  6       FALSE     1     6          0         35
## 7   6.63149501  5.36850499  7       FALSE     1     7          0         36
## 8   7.53505551  6.46494449  8       FALSE     1     8          0         36
## 9   8.15720719  5.84279281  9       FALSE     1     9          0         23
## 10  7.33832123  6.66167877 10       FALSE     1    10          0         36
## 11  7.42440009  6.57559991 11       FALSE     1    11          0         43
## 12  8.59464084  7.40535916 12       FALSE     1    12          0         32
## 13 10.30218091  9.69781909 13       FALSE     1    13          0         43
## 14 11.98835832 10.01164168 14       FALSE     1    14          0         38
## 15 14.34135852 13.65864148 15       FALSE     1    15          0         48
##      xmin   xmax xid newx new_width weight colour size alpha shape linetype
## 1   0.625  1.375   1    1      0.75      1 grey20  0.5    NA    19    solid
## 2   1.625  2.375   2    2      0.75      1 grey20  0.5    NA    19    solid
## 3   2.625  3.375   3    3      0.75      1 grey20  0.5    NA    19    solid
## 4   3.625  4.375   4    4      0.75      1 grey20  0.5    NA    19    solid
## 5   4.625  5.375   5    5      0.75      1 grey20  0.5    NA    19    solid
## 6   5.625  6.375   6    6      0.75      1 grey20  0.5    NA    19    solid
## 7   6.625  7.375   7    7      0.75      1 grey20  0.5    NA    19    solid
## 8   7.625  8.375   8    8      0.75      1 grey20  0.5    NA    19    solid
## 9   8.625  9.375   9    9      0.75      1 grey20  0.5    NA    19    solid
## 10  9.625 10.375  10   10      0.75      1 grey20  0.5    NA    19    solid
## 11 10.625 11.375  11   11      0.75      1 grey20  0.5    NA    19    solid
## 12 11.625 12.375  12   12      0.75      1 grey20  0.5    NA    19    solid
## 13 12.625 13.375  13   13      0.75      1 grey20  0.5    NA    19    solid
## 14 13.625 14.375  14   14      0.75      1 grey20  0.5    NA    19    solid
## 15 14.625 15.375  15   15      0.75      1 grey20  0.5    NA    19    solid
# Pairwise comparison of all sites
stat.test <- cont_snotel_site_yr_apr1 %>%
  mutate(US_L3CODE = as.numeric(US_L3CODE)) %>% 
  arrange(US_L3CODE) %>% 
  pairwise_t_test(newSWE_days ~ US_L3CODE, p.adjust.method = "bonferroni") %>% 
  add_significance(cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
                   symbols = c("****", "***", "**", "*", "ns"))
  
level1 = c(4,5,9,11,13,15,16,17,18,19,21,23,41,77,80)
level2 = c(5,9,11,13,15,16,17,18,19,21,23,41,77,80)

pairwise_plot <- ggplot(stat.test, aes(x = factor(group1, levels = level1),
                                       y = factor(group2, levels = level2),
                                       fill = p.adj.signif))+
  geom_tile(col = "black") +
  scale_fill_brewer(labels = c("1e-04", "0.001", "0.01", "0.05", "ns"), 
                    palette = "Reds", direction = -1)+
  labs(y = "EPA Level III Ecoregion", x = "EPA Level III Ecoregion", fill = "Sig. Level") +
  theme_bw()

pairwise_plot

ggplotly(pairwise_plot)
#test <- ggplot_build(pairwise_plot)$data[[1]]

# boxplot and pairwise 
plots <- plot_grid(
  spring_snow_days + theme(legend.position="none"),
  pairwise_plot,
  align = 'vh',
  labels = c("a.", "b."),
  hjust = -0.4,
  nrow = 1)
## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.
plots

legend <- get_legend(spring_snow_days + guides(fill = guide_legend(nrow = 5)) +
               theme(legend.direction = "horizontal",
                     legend.justification="center",
                     legend.box.just = "bottom"))

plot_grid(plots, legend, ncol = 1, rel_heights = c(1, .3))

# Elevation vs new SWE count (March 1)
elev_newsnow <- ggplot(cont_snotel_site_apr1, aes(x = med_spring_days, y= elev, 
                                                  color = reorder(eco_sntl,median_spring_SWE,na.rm = TRUE))) +
  geom_point() +
  #geom_smooth(method = "lm")+
  #stat_smooth(method = "lm", col = "black")+
  labs(x = "Median Days of Increased SWE", y = "Elevation (m)", color = "") +
  scale_x_continuous(breaks = seq(0,30,10), limits = c(0,30), expand = c(0.01, 0.01)) +
  scale_y_continuous(breaks = seq(500,4000,1000), limits = c(500,3750), expand = c(0.01, 0.01)) +
  theme_bw()
elev_newsnow

ggplotly(elev_newsnow)
#ggplot_build(elev_newsnow)$data

plot <- plot_grid(
  spring_snow_days + theme(legend.position="none"),
  elev_newsnow + theme(legend.position="none"),
  align = 'vh',
  labels = c("a.", "b."),
  hjust = -1,
  nrow = 1)

plot

legend <- get_legend(spring_snow_days + guides(fill = guide_legend(nrow = 5)) + 
               theme(legend.direction = "horizontal",
                     legend.justification="center",
                     legend.box.just = "bottom"))

plot_grid(plot, legend, ncol = 1, rel_heights = c(1, .3))

#Determine the number of SWE events above a temperature threshold
summary_stats <- cont_snotel_spring_apr1 %>%
  #filter(site_id == "1000" & (date >= "2001-04-01" & date <= "2001-04-15")) %>% 
  summarise(n_above_mean = sum(newSWE > 0 & temperature_mean > 0, na.rm = TRUE), 
            n_above_min = sum(newSWE > 0 & temperature_min > 0, na.rm = TRUE),
            n_below_mean = sum(newSWE > 0 & temperature_mean <= 0, na.rm = TRUE), 
            n_below_min = sum(newSWE > 0 & temperature_min <=0, na.rm = TRUE),
            n_newSWE = sum(newSWE > 0, na.rm = TRUE),
            n = n())

above_min_pct = summary_stats$n_above_min/(summary_stats$n_newSWE)